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ABSTRACT 


An acoustic tomography array consisting of six transceiver moorings was 
jointly deployed by Woods Hole Oceanographic Institution and Scripps 
Institution of Oceanography in the Greenland Sea during the second half of 1988. 
Two of the primary objectives of this thesis are: (1) to set up and test a stochastic 
3-D inversion code for the Greenland Sea Acoustic Tomography data analysis; 
and (2) to evaluate the performance of the acoustic system through resolution 
and variance analyses. In acoustic tomography, the sound speed perturbation 
field is estimated from measured acoustic travel time perturbation data. A unique 
sound speed perturbation estimate can be constructed using the Guass-Markoff 
theorem. However, the theorem requires the specification of the covariance of 
the sound speed perturbation field, which is generally not exactly known. Via 
computer simulation, we examined the sensitivity of the estimate to uncertainty 
in the sound speed field correlation specified. In addition, we also examined the 
effects of an increased random experimental noise level and a change in array 
geometry due to mooring failure on the estimate. The three major results are 
that: (1) the estimate is less sensitive to a positive uncertainty in correlation 
length than to a negative uncertainty in an ocean volume containing large 
structures, while it is more sensitive to a positive uncertainty than to a negative 
uncertainty in an ocean volume containing small structures; (2) the estimate 
error is primarily bias error rather than random error; and (3) the failure of a 
mooring causes a large increase in RMS error in regions no longer containing 
acoustic rays, but it results in an increase in RMS error of only 25% in regions 


which still contain acoustic rays. 
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I. INTRODUCTION 


A. OCEAN ACOUSTIC TOMOGRAPHY 

Ocean Acoustic Tomography is a method used to monitor the mesoscale 
ocean variability (which is the oceanic analog of atmospheric weather) and it was 
introduced by Munk and Wunsch (1979). This technique is analogous to the 
medical X-ray procedure known as Computer Assisted Tomography (CAT) 
(Figure 1-1.a). Roughly speaking, tomography exploits the fact that the ocean is 


"transparent" to acoustic rays to remotely sense the properties of an ocean 


region. 
X-ray source ocean front 
y Onc ae a a 
- 7 Se 
body - 7 i — 





acoustic @=---- 
transceiver 





(a) Medical CAT (b) Ocean Acoustic Tomography 


Figure 1-1: The Comparison of Medical CAT and Ocean Acoustic 
Tomography. 


In practice, a number of acoustic transceivers are deployed at positions chosen to 
allow for coverage of an ocean volume of interest (such as a region containing 


mesoscale eddies or a frontal system) (Figure 1-1.b). The most common 


application of tomography is for estimating the perturbation of the sound speed 
field from a set of measured acoustic travel time perturbations. The 
perturbations in sound speed are assumed to be so small that the perturbations in 
acoustic travel time between each pair of transceivers are linearly related to the 
sound speed perturbations. The modeling of the travel time perturbations due to 
the sound speed perturbations is known as the forward problem. Once the 
forward problem is solved, inverse methods which are widely used in 
geophysical research (Backus and Gilbert, 1967) are applied to the travel time 
data for the reconstruction of the the sound speed perturbation field. 

Ocean Acoustic Tomography offers several advantages over conventional 
hydrographic surveying method. These advantages are pointed out by Chiu 
(1978): (1) the system can be implanted in the ocean on a semipermanent basis to 
allow for continuous observation; (2) it is not affected greatly by weather 
conditions; (3) it has high temporal resolution; (4) it can cover an extensive 
volume of the ocean interior and probe the different parts simultaneously; and 
(5) only a few moorings are needed, thus minimizing the effort in deployment 
and maintenance. 

Since the first successful experiment (the 1981 Three-dimensional Mesoscale 
Experiment), additional tomography projects have provided measurements of 
mesoscale eddies (Cornuelle et al, 1983), planetary waves (Chiu et al, 1987), 
currents (DeFerrari et al, 1986), internal waves (Stoughton et al, 1986), basin 
mode oscillations (Bushong, 1987), and surface waves (Lynch et al, 1987). In the 
future, monitoring of large-scale ocean dynamics on a global basis may be 


achieved using cross-basin transmissions. 


B. GREENLAND SEA PROJECT OCEAN ACOUSTIC 
TOMOGRAPHY 


The Greenland Sea Project (GSP) is a plan developed by the international 
Greenland Sea Science Planning Group which was appointed by the Arctic Ocean 
Sciences Board (AOSB). The overall goal of this five-year program (from 1987 
to 1992) defined by AOSB is to understand the large scale, long-term 
interactions among the air, sea, and ice in the Greenland Sea. The primary 
region of the study is bounded by Fram Strait to the north, Spitsbergen and the 
Mohn Rise to the east, the Greenland-Jan Mayen Ridge to the south, and 
Greenland to the west (Greenland Sea Science Planning Group, 1986, pp. 1-7). 

The plan is designed to study the following ocean dynamics: (1) the seasonal 
and interannual variability of the sea ice cover; (2) ocean ventilation and 
convection of the deep water; (3) ocean circulation and mixing; (4) atmosphere 
energetics; and (5) biological processes. The Ocean Acoustic Tomography Array 
is a component used in GSP to monitor the process of ocean ventilation and 
convection in the Greenland sea central gyre (Greenland Sea Science Planning 
Group, 1986, pp. 1-7). It is the process of ventilation and convection that gives 
the Greenland Sea central gyre the ability to affect the oceans throughout the 
world. 

An acoustic tomography array consisting of six transceiver moorings with a 
pentagonal geometry was jointly deployed by Woods Hole Oceanographic 
Institution and Scripps Institution of Oceanography in the Greenland Sea during 
the second half of 1988. Figure 1-2 shows the tomography array position and the 
GSP area. 
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Figure 1-2: Acoustic Tomography Mooring Array Position and Geometry 
Configuration. 


For convenience of calculation, the given geodetic position (latitude and 
longitude) of each array element has been translated into an xy position using a 
Transverse Mercator (TM) projection. To do this, the longitude of the central 
mooring (array #6), was chosen as the central meridian of the TM projection. 
The false origin was set at 73°N and shifted to 120 km west of the central 
meridian. This coordinate system gives only slight position distortion at the edge 
of the array. The ray paths can then be calculated in the xy planar coordinates 
rather than in geodetic coordinates. TABLE 1-1 shows the coordinate conversion 
for the moorings as well as the depths of the acoustic sources and vertical 


receiver arrays. 


TABLE 1-1: COORDINATE CONVERSION OF THE ACOUSTIC 
MOORINGS (USING WGS72 SPHEROID) 


Mooring Lat. Lon. Xx y S-depth* R-depth™ 
(meter) (meter) (meter) (meter) 
TS4AT.IN  1°04.7'W 150726 312598 94.5 145.3 
2 T453.3N 1°24.8'E 225124 213940 D504, 145.4 
3 74°00.4'N — 1°05.2'W 154249 112659 94.6 145.1 
4 T4°18.9N 4°59.0'W 36008 148715 94.6 117.0 
» feo. 2 No) 15.9 W 33882 272300 94.8 Nee 
6 T4’°54.0N 2°12.0'W 120000 212040 oe 145.7 


* The depth of acoustic source. 


** The depth of vertical acoustic receiver array. 


C. THESIS OBJECTIVE 

The purpose of this research is to develop computational tools in preparation 
for analysis of the tomography data from the GSP. Specifically, our objectives 
are: 

¢ to set up and test a stochastic 3-D inversion code for the Greenland Sea 
Acoustic Tomography data analysis; and 

¢ to investigate the sensitivity of the sound speed perturbation estimate to 
Our uncertainty in the sound speed field correlation, changes in the random 
experimental noise level, and changes in array geometry due to mooring failure. 
The second objective is accomplished via computer simulation. 

In Chapter II, we discuss the behavior of an acoustic ray in an 
inhomogeneous moving medium. A linear sound speed profile is taken as the 
reference state of the sound speed field. Based on this reference state, a 
numerical 4th order Runge-Kutta integration method was used to calculate the 
paths of the eigenrays for establishing the forward acoustic model. We then 
discretize the forward model by dividing the ocean volume into 500 boxes in 
order to cast the problem into a matrix form for the computer simulation. 

In Chapter III we present a three-dimensional stochastic inverse method 
which is the distribution-free Gauss-Markoff estimator. Due to insufficient 
experiment data, the inverse problem is underdetermined. In this stochastic 
approach, a priori information is specified in the covariance matrix of sound 
speed perturbations. The covariance gives additional constraints to the system 
and therefore a unique Solution is obtained. Two measures, RMS error and 
resolution length, are used to quantify the performance of the estimator at each 


box location. 


The ocean is a dynamic and inhomogeneous environment. It is difficult to 
incorporate an exact sound speed perturbation covariance matrix as a constraint 
since we do not have enough statistical information at this tme. Therefore, an 
approximate sound speed perturbation covariance matrix is generally used. When 
the covariance is inexact, the estimator 1s suboptimal. In Chapter IV, we vary the 
assumed correlation length (which is used to construct the sound speed 
perturbation covariance matrix for the estimator) to determine the sensitivity of 
the system to uncertainty in correlation length. We also study the effects of 
experimental random noise and failure of array elements on the estimate on the 
estimate. 

In Chapter V, we present a summary of our research, along with results and 
conclusions. Furthermore, we propose a criterion for designing estimators for 
the analysis of the GSP data. Finally, we make recommendations for improving 


our research. 


Il. ACOUSTIC FORWARD MODELING 


A. TRAVEL TIME 

In acoustic tomography, the sound speed and flow fields are reconstructed 
from travel time measurements. The corresponding forward problem is, 
therefore, to find eigenrays, 1.e., those rays which emitted by the source that are 
intercepted by the receiver, and to establish the relation between the 
measurement and the unknown fields. Each eigenray has a unique launch angle, 
and a unique path through the ocean, thereby sampling the sound speed field and 
flow field, at different locations (Cornuelle, 1983, pp. 41). 

The geometric approximate travel time along a ray path, T(t), can be 
calculated by integrating the ray slowness along the path. In the presence of a 


current v(r,t), T(t) is given by 


T(t)= ds 


nan? 


en c(r,t) + v(r,t)°T 


(2a 

where c(r,t) is the sound speed at time t and position r, v(r,t) is the ocean 
eA ; 

current velocity, T is a unit vector tangent to the ray, and s is the arc length 

along the path. The travel time is changed by the sound speed perturbation field, 

dc(r,t), which is the deviation from the reference sound speed, c,(r). The 


reference sound speed can be the overall space-time average sound speed. Thus 


the travel times in a reciprocal transmission can be expressed as 


T'=T,+8T = as - 
aT c,(r,t) + Sc(r,t) + v(r,t) *t (2.2) 
and 
b 
T’ =T,+8T = eS 
a c.(r,t) + Sc(r,t) — v(r,t) *t (2.3) 


where the superscripts f and b refer to forward and backward transmissions, 
respectively, and ST is the perturbation of the reference travel time T,. 
: ; —2 
In most ocean environments, dc/c, is on the order of 10 ~. Thus we can 


approximately linearize the reciprocal travel time perturbations as 


f ar os 
qa | Sect) + v(r,t) : fe 





co(r) 
path (2.4) 
and 
a 7 Sc(r,t) ; v(r,t)*T ms 
C,(r) 
path (5) 
Taking the sum of Eq. (2.4) and Eq. (2.5), one obtains: 
f b ; 
c= oT +0T =~— sb ds, 
: co(r) 
path (2.6) 


while taking the difference of Eqs. (2.4) and (2.5), one obtains: 


o~ 


f b ary 
~ oT -d8T v(r,t) °T 
= Res ds 


2 co(r) 
path (25 


where 5T’, half of the summation of the forward and backward travel time 
perturbations, is linearly related to the sound speed, and OT , half of the 
difference of the forward and backward travel time perturbations, 1s linearly 
related to the current. In this thesis our focus 1s on the estimation of the sound 
speed perturbations only, and we won't be dealing with Eq. (2.7) at all. 

In order to express the travel time in a vector form the continuous integral 
in Eq. (2.6) was discretized by dividing the ocean into small boxes, in which the 
sound speed perturbation was assumed constant. Figure 2-1 shows the horizontal 
transceiver array configuration and the corresponding box geometry in a 
horizontal slice. 

In this project we discretized the ocean volume by 500 boxes (10x10 squares 
horizontally and 5 layers vertically); the limitation of computer memory space 
in the microVAX dictated this decision. 

After the discretization of the sound speed perturbation field into a vector 
oc, the vector OT containing all the eigenray travel times can be expressed ina 


matrix-vector form as 


OT = A oc, (2.8) 


where 


(2.9) 


10 


is the element in the ith row and jth column of the matrix A and is equal to the 


integral of als: along a segment of the /th acoustic ray path in the jth ocean 


box. 





Figure 2.1: The Horizontal Transceiver Array Configuration and the Box 
Geometry System Used. 


B ACOUSTIC RAYS IN| AN INHOMOGENEOUS MOVING 
MEDIUM 


A numerical ray-trace algorithm for the acoustic rays in an inhomogeneous 
moving medium was used in our study. The procedure was developed by Chiu 


(1985). A summary of the theory will be given in the following discussion. 


1. Ray Path 
It is well known, in acoustic ray theory, that the acoustic rays in a 
motionless medium with space dependent index of refraction, n(r), is described 
by the eikonal equation: 
| Veo(r) , =n(r), (2.10) 
where ®(r) is the eikonal function (1.e., acoustic phase) defining the wave fronts, 
r = (x,y,z) is the position vector, and V@(r) points to the direction of 
propagation of acoustic field energy. However, when the inhomogeneous 
medium is moving with a velocity given by v, the governing eikonal equation for 
acoustic wave fronts becomes (Ugincius, 1970) 
2 


lver)| = = n(r) z6 —yv- ), 


(2.11) 
Ugincius (1970) has derived a second-order vector differential equation 
governing the ray paths based on the eikonal equation Eq. (2.11). The equation 


for ray paths is 
& (Nr') = (r"*V)V + r'x(VxV) = VN, 


(2.12) 
where 
Nas eee (2.13) 
S=y l-v +(r's vy)’, (2.14) 
and 
7 n(S —r'ev) 
(1-o'}s_ (215) 


| 


and where N and V are functions of both the position r of a point as well as the 
ray direction r' through that point, and v =v/c is nondimensional medium 
velocity and it has a magnitude equal to v. Note that, we have used ‘ and " to 
denote the first and second derivatives with respect to s, respectively. 

In the ocean, the current velocity is small compared to the speed of 
sound. Therefore, Eq. (2.12) can be simplified by eliminating the second order 


terms involving v in N and V so that 


N=n , Ve=n(l—-r'ev)v, (2.16) 


and 


(r"eV')V = 0. eal) 


By doing so Eq. (2.12) reduces to 


d ? ' ' — 
= (ne') + 4°x[Vxn(1 - r'sv)v |= Vn. 


(2.18) 
By replacing r' with dr/ds in Eq. (2.18) we have 
ddr. dr dr 
qs (gy) + $bq Van(1 SE + vy |= Vo. (2.19) 


In this thesis we neglected the ray curvature in the horizontal plane and 
let the reference sound speed profile to be a function of depth z only. This 
restriction implies that the reference ray paths are confined to lie on vertical 
slices normal to the xy plane (Ziomek, 1985, pp. 232). In Figure 2-2 a ray path 
with such a restriction is shown. On a vertical slice (i.e., a Range-Depth Rz 


plane), the equations governing the planar rays on this plane are 
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Fa (ne = gsr) Ole 
OR (2.20) 
and 


qe dz a a ona 
i OG ~ Ger)? = 3, 





(2.21) 


A oo" ; 
where R and 7 are the unit vector in R and z direction, respectively. 






z (depth) 


Figure 2-2: The Ray Path Restricted in a Vertical Slice Normal to xy 
Plane. 


Since the reference sound speed profile is taken to be a function of 


depth only, we can simplify the differential equation (Eq. (2.20)) to get 


dR cos8,— v.+n(z)vp 
ds = nz = cos(8) 
(2.22) 


2 
=+ 4 / | ——%@ —__| - 1 =tan@). 
cosO . — Uo + n(z)UR (2.23) 


where v, is the medium speed at position (R,, Z,) in the ray path direction. As 


OT 


|e 


shown in Figure 2-2 0, and 6 are the initial launch angle at position (R,, z,) and 
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ray angle at position (R, z), respectively; Eq. (2.23) can be integrated in either 
the R or z directions to get the ray paths. In our raytracing we divided each 
range connecting each of the transceiver pairs in the Greenland Sea into 1,000 
steps and integrated Eq. (2.23) in the R direction for the corresponding depth 
values. Integrating in range gives, for each ray path, depth of the trajectory as a 
function of range. The function has a one-to-one correspondence. 
2. Numerical 4th Order Runge-Kutta Integration Method 

An accurate ray path can be calculated using a well known numerical 

4th order Runge-Kutta integration method. From Eq. (2.23). the integral for 


which we need to compute Is: 


R 


2 
pene ny) |_| - 1 ar. 
cosO,— U9 + N(z)UR 


Ro (2.24) 
The numenical 4th order Runge-Kutta integration method is given by 


the following formulae (Gerald, 1989, pp. 358): 


R 


Znet=Znt f f{R, Z(R)| GR = 2_ +2 (ky+ 2k>+ 2k3+ky), 


Ra (si 

where 

2 
f (R, z(R)] =+ ca NS ee ee 
cos8 — UV 5tN(z)U R (2.25.2) 
k, =hf(R,, z,), (2.25.3) 
kp =hf(R,+th, z+ 5k), (2.25.4) 
k3=h f(R,+ 4h, z+ 7k), (2.25.5) 
and 


lis 


k, =hf(R, +h, z, +k): (2.25.6) 

A step size h has been selected to limit the numerical error to a 

tolerable size. The global error (1.e., discretization error) accumulated along the 

entire interval R by the 4th order Runge-Kutta method is (Gerald, 1989, pp. 
358) 


O(h’) = Rh’ f{E.z6)] . O<E<h. 
180 (226) 


The error cannot be exactly determined because the position (R, z) = [& , z( &)] is 
an unknown with € bounded by the interval (0, h]. A standard way to determine 
whether the z values are sufficiently accurate is to compare the value computed 
using a step size of h with the value calculated using the half of h. If this gives 
only a change of negligible magnitude, the results are accepted; if not, the step 1s 
halved again until the results are satisfactory (Gerald, 1989, pp. 358). 

3. Turning Point 

AS an acoustic ray travels through an ocean volume it will be refracted 
upward or downward. When the ray angle goes to zero, the ray will start to 
bend up or down. The position of this zero ray angle point is called the turning 
point. 

When a ray is traced to a position which is less than a step size away 
from a turning point, we use a linear gradient approximation to calculate the 
depth of the ray at the terminal of that step. This approximation near a turning 
point is needed because the function for which we integrate will no longer be the 
same after the turning point. There will be a sign change in the function. The 
curvature K can be used to calculate the ray path around turning point in the 
linear gradient case. Ugincius (1970) has derived the following equation of 


Curvature: 
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a es ee ey gata 1) 
—¢ dz ds dz ds QT) 


Let the gradient of sound speed dc/dz be g, and the gradient of medium speed 
dv/dz be g,, Eq. (2-27) becomes 


l 
K =—ig,(v —cos®) + g,(2cos8 v — 1)}. 
= |e ) + gol ) ae 


Since g, and v are very small, the term involving the product of g, 


and v is negligible. Thus Eq. (2.28) can be approximated by 


Igo COST) = 9»)| 
K = - 


C (2.29) 


The radius of curvature ® 1s the reciprocal of k, 1.e., 


R= 5 en SS 

lado = ase) — 9») (2.30) 
Under the assumption of constant gradients near a turning point, the 
local radius ® is also a constant, which result in a circular ray path, locally. If 
the radius ® were negative, the radius would curve upward and vice versa 
(Kinsler, 1982, pp. 401-402). In the following discussion R will be referred to as 

the magnitude of the radius of curvature. 
The geometries of segments of ray paths through turning: points are 
shown in Figure 2-3 and Figure 2-4, where @, is the ray angle before the turning 
point and 9, is the ray angle after the turning point. The corresponding depth 


increment can be equated as 


Ge 0;— CX COS0 5 COSU,,), (2.31) 


17 


~ 
step size h = 





z (depth) 


Figure 2-3: The Geometry of a Downward Turning Ray Path. 
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Figure 2-4: The Geometry of an Upward Turning Ray Path. 
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where 
{2 erate 
Rcos8, = i (h = R sind » / (2232) 


A] ee 


and "+" is for downward curve while "—" is for upward curve. From Eq. (2.31) 


the depth of the turning point is simply given by 


Zw = Z;+ R(1—cos8,). (2.33) 
4. Surface Reflection 

Due to the upward refracting nature of the Greenland Sea sound 
channel, surface reflection of rays is of importance. The ray path is 
discontinuous at the point of reflection and the equation for rays Eq. (2.24) is not 
valid here. We need to develop a procedure to calculate the local path trajectory 
near a reflection point in an other way. The same linear gradient assumption, as 

used for the calculation near a turning point is applicable here again. 
Similar to Eq. (2.31), the depth increment in a range step within which 


a surface reflection occurs 1S 


Za 2, — = (XK COSO> eR COSO ,). (2.34) 


In the case of convex path segment (figure 2-5), Rcos®, is given by 


/ 2 
Rcos0, = R - (2Rsin8 — RsinO, —h) (2.35) 
= D 
RsinO = “V R —(Rcos8, + Z)) . 


In the case of a concave path segment (Figure 2-6), Rcos®, 1s given by 


and 


(2.36) 
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Figure 2-5: The Geometry of a Surface Reflected Convex Path Segment. 
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Figure 2-6: The Geometry of a Surface Reflected Concave Path Segment. 
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2 
Rcos0, = V R — (2Rsin® — KsinO, +h) Og 
a 2 : 
Ksin@= VY KR —(Keos8,—Z;) . (2.38) 


Having all the equations coded in FORTRAN, predictions of ray paths as well as 


and 


eigenray finding were possible. 
5. Ejigenrays Finding 
Since the sound-speed profile in the central Greenland Sea gyre is very 
nearly adiabatic below the surface layers, a linear and range-independent sound 
speed profile was selected to represent the reference state c,(z), as shown in 
Figure 2-7 (left). We then use this reference state to calculate the ray pattern of 


transmission in the acoustic forward problem. 
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Figure 2-7: Typical Sound Speed Profile and Ray Path with The Ray Angle 
from -15° to 0° in The Greenland Sea Project Area. 


In Figure 2-7 (right) the ray paths shown are the typical propagation 
pattern in the arctic environment. The rays with fewer loops, that reach the 
deeper layers, are the faster paths for the acoustic energy. 

In order to search for rays that reach the receiver (1.e., eigenrays) we 
need to shoot rays with a range of launch angles. Given an angular interval 
within which all possible eigenrays lie, we shall be able to determine the 
eigenrays by looking for the intersections between the arrival depth curve and 
the receiver depth line as illustrated in Figure 2-8. Three eigenray patterns 
associated with 3 different ranges are shown in Figure 2-9. We only choose 6 to 
7 out of as many as 30 eigenrays that sample the ocean from the surface to 1 km 
depth. Figure 2-10 shows the selected eigenrays for the three ranges. There are a 


total of 91 eigenrays will be used to conduct our sensitivity study. 
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Figure 2-8: Arrival Depth Curve for Eigenrays Finding. 
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Figure 2-9: Typical Eigenray Pattern for Three Different Range as in (a) 
Range = 123.6 km, (b) Range = 105.2 km, and (c) Range = 200.1 km. 
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Figure 2-10: Typical Selected Eigenray Pattern for Three Different Range 
as in (a) Range = 123.6 km, (b) Range = 105.2 km, and (c) 
Range = 200.1 km. 
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Ill. STOCHASTIC INVERSE METHOD 


In the forward problem the ray travel times are modeled using Eq. (2.9). In 
this chapter, we discuss the inversion of Eq. (2.9) using a Gauss-Markoff 
estimator to obtain an optimal estimate of the three dimensional sound speed field 
in the ocean volume monitored by the acoustic array. The estimation error and 
system resolution will be analyzed in the next chapter in an effort to quantify the 


performance of the array. 


A. ESTIMATION OF SOUND SPEED FIELDS. 
1. The Gauss-Markoff Estimator 
In the Gauss-Markoff stochastic inverse method, both the data 6T and 
the unknown field 6c are assumed to be random vectors. The forward model 


relating the data and the unknown field is 


oT=AX+e, (3.1) 
where e is the experimental noise which corrupts the travel time measurement; e 
is assumed to be uncorrelated with dT and dc. 

The system Eq. (3.1) is highly underdetermined, and thus without 
additional constraints other than just the data, the system admits infinite number 
of solutions for 6c. In stochastic inverse methods a unique optimal linear 
estimate of the unknown parameters, Oc, can be constructed by incorporating a 
priori knowledge of the parameters in a covariance matrix. 

The well known Gauss-Markoff estimator is chosen here because it 


requires no knowledge of probability densities. This so called distribution free 
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property is the most important feature of the Gauss-Markoff estimator (Liebelt, 
1967, pp. 136). The estimate Sc Satisfies a minimum mean square error 


criterion 


Cem ) = minimum, 


e2) 
Following Liebelt (1967) and Chiu (1987), the Gauss-Markoff estimate 
is given by 
Se= CAC,’ ST, (3.3) 
where 


a | 
=| T~-1 
Galea C al (3.4) 


. A ; 
is the covariance matrix of the error € = Oc— Oc in the estimate, and C, and C5, 
are the covariance matrices of noise e and the unknown parameters &, 


respectively. The construction of the estimate requires finding the inverse of the 
' -1 Te 
matrix Cs. + A C.A. 


The trading of system resolution for stability in the optimal estimate can 
1 


be revealed by a singular value decomposition of the matrix C, zs ond such that 
| : : 
CAA oH = 
: Se (Gi) 


where the diagonal elements A; of the matrix A are the associated singular values, 
and the columns u, and v; of U and V are the left and right singular vectors, 
respectively. 

The matrix, Eq. (3.5), is the operator associated with a 


nondimensionalized version of the forward model: 


924 


1 


2S = =i 
C. ar-|c, acy C5 + C. c 


(3.6) 
which is transformed from Eq. (3.1). Using Eq. (3.5), Chiu et al (1987) have 


shown that the minimum mean square error estimate can be expressed as 


1 





2\ a 
I+A]/ ALU C oT 





say 
C50 (3.7) 
or equivalently as 
als k Xr. ee 
Ghee = »y : luc? ar Vi 
i=l 
eee | (3.8) 


The vectors v; are the base vectors occupying the solution domain such 
that any solution can be expressed by a weighted sum of these vectors. A singular 
value A; much small than one is associated with a singular vector v; that models 
a highly unstable and oscillatory function. The linear estimator downweights 


these oscillatory function to stabilize the estimate. 


B. ERROR AND RESOLUTION 
In the previous section, we introduced an estimator which solves the inverse 
problem. In this section, we derive measures which are used to quantify system 
performance. 
1. Error of the Estimate 
The error covariance matrix of the estimate was expressed in Eq. (3.4), 


which can also be equated, using Eq. (3.5), as 


1 


a (3.9) 


1 


C,=C? 


fea) 
ge LE- VA + A KV Ee 





The diagonal terms of this matrix are the mean square errors of the estimate at 
each box. 

The error of the estimate, Soe dc, has two components, bias and 
random error. The bias, b = (Sc) — 6c, results from an insufficient number of 
data, and the random error, A(Se ) = Sc _ (8c) is caused by the random 


experimental noise. Since these two components are statistically independent, C, 


can be expressed as 


1 

C.= (bb ) + Cage (3.10) 
where 

7 jee 

Caisey=CeA C, AC, (3.11) 

or equivalently, 
=] =] 
5c bast hesT vy 
Case) = C5. VII+A A \I+A V Ce. (3.12) 


An expression for (bb' ) can be obtained by letting the random error 


covariance C,/s,) approaches zero, i.e. by letting all the eigenvalues 4; approach 


infinity in Eq. (3.9). The result is 


— 


Te ; ie ; 
bb = I~-VV | C.. 
(ob") =c% [t-vv"| ch 2 
2. Resolution 
We define a symmetric, nxn matrix 
=) 
2 T 


as the resolution matrix. Using the definition of R, Eq. (3.9) becomes 
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C.= Ci LI-R] Co. ag 

If the resolution matrix 1s an identity matrix (that is, R=I), then the 

efror 1s zero. On the other hand, if the resolution matrix is not an identity 
matrix, for example, if it has nonzero diagonal elements then the error 1s 
nonzero. This implies that the mean square error is intimately related to system 
resolution. The ith row of the resolution matrix, (r.}°, is defined as the 
resolution kernel of the ith box and it describes how much neighboring boxes 
contribute to the error of the estimate in the ith box, or how well the field at the 
ith box can be resolved (Menke, 1984, pp. 61-68). Figure 3-1 illustrates the 


structure of a typical resolution matrix. 


Figure 3-1: Plot of Selected Rows and Columns of The Resolution Matrix 
R 


If the ith resolution kernel has a single spike centered at the diagonal, 


then the ith box is perfectly resolved. If the peak 1s very broad, the ith box is 
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poorly resolved and has a large error due to the fact that the estimate at the ith 
box is an average of the neighboring field. In an other word the energy in the 
estimate is spread to neighboring boxes due to the poor resolving power of the 
system. 
Two measures of the resolution of the estimator have been proposed. 
One is the so called resolution "spread" (Miller, 1989, pp. 333), which measures 
the difference between the resolution matrix and an identity matrix via the 
expression 
Spread = y y OS | mF 
a we (3.16) 
The second measure is the so called minimum resolution length (Chiu, 
1987), which measures the local resolution at each box location. The resolution 
length is defined to be the distance at which the resolution energy falls to half of 
its peak value. To be more precise the minimum resolution length at the ith box 
is defined as the square root of the second central moment of energy distribution 
in the ith resolution kernel ({r,}'). The minimum resolution lengths in the three 


spatial directions may be expressed as 





aan) > Y Siew a) * ni Giydie | 









IX Jy Jz ni (3.17.a) 
Hy(ix,ly,1Z) = SY ¥ tGy- ly) dy ] Oy 
aes Yael ri (3.17.6) 


and 
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Uixiinsal > SS Eee * HGrdyia) | 





x Jy ye ; (3:1 ie) 
where E, is the total energy of the ith resolution kernel r; . 
=> d axiyia). 
yd ae (3.17.d) 


In Eq. (3.17.a-b) #4, and H, are the minimum horizontal resolution lengths in the 
x and y directions, respectively, while V is the minimum vertical resolution 
length. Physically, the minimum resolution lengths determine the minimum eddy 
size that can be resolved adequately by the monitoring system. Note that we have 
expanded the row index / and column index / into three dimensional box indices, 
(1x, 1y, 1Z) and (x, Jy, jz), in Eqs. (3.17.a-d). The expansion was done according 


to the following equations. 


i = (ix-1)Xnyxnz + (iy-1)xnz + iz = 1, 2, ...... n, (3.18.a) 
= (jx-1)xnyxnz + (jy-1)xnz + jz = 1, 2, ...... n (3.18.b) 

with 
n = nxxXnyxnz, (3.18.c) 
where 1X; {xX —) lee nx, is box indices in the x direction; iy, jy = 1, 2, 
oe ny, is box indices in the y direction; iz, jz = 1, 2, ......nz, 1s box indices in the 


z direction. With nx = 10, ny = 10, and nz = 5, the total number of boxes is n = 
500. 
From the above discussion in this chapter, It has been found that the 


RMS error and resolution length do not depend on the experimental data OT, but 


only on the sound speed perturbation covariance matrix C;., travel time error 


jy 


covariance matrix C,, and the transfer function A (which depends on the array 


geometry and the number of rays). 
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IV. RESULTS OF SENSITIVITY STUDY 


In the last chapter, we outlined two ways to evaluate the performance of the 
Greenland Sea tomography array by examining the statistics of the error in the 
estimate and the resolution of the system. As mentioned the RMS error of the 
estimate and the system resolution measures depend on the covariance matrix of 
sound speed perturbation C;., the covariance matrix of noise C,, the array 
geometry, the number of eigenrays but not the data OT. Therefore, even without 
using actual or synthetic measurements, the performance of the array and ray 
path geometry can be evaluated given the covariances. By varying the ocean 
correlation length (which determines C;.), the rms value of noise (which 
determines C,.), and the array geometry (which determines A), we have 
examined the changes in system performance as these ocean and acoustic 
parameter vary. The results are discussed in Section IV-B. 

By incorporating statistical information concerning the covariance of the 
field, the indeterminacy of the unknown field, dc, is eliminated. Since the 
covariance matrix Cs. is the a priori information that we supply to the 
estimator, we are particularly interested in determining the sensitivity of the 
estimate to the uncertainties in the correlation lengths of the sound speed 
perturbation field. This sensitivity study was accomplished through inversions of 
synthetic data generated in the computer. The results are discussed in Section IV- 
©: 

In the following section, we first discuss the method we use to simulate sound 


speed perturbation fields and travel time data in the computer. 
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A. COMPUTER SIMULATION OF MESOSCALE SOUND SPEED 
FIELDS AND TRAVEL TIME DATA 


The temperature, sigma-t, and sound speed profiles obtain by a CTD 
(Conductivity, Temperature, and Density) cast in the Greenland Sea (Worcester 
and Howe, 1989) is shown in Figure 4-1. We have superimposed on this data a 
linear profile (dash line) which we have chosen as the reference sound speed 
profile. In the simulation work, we take that the perturbations of sound speed 
occur only in the water column shallower than 1,000 m. This is not a bad 
assumption as indicated by CTD profile obtained by Worcester and Howe. For 
simplicity, we also neglect ocean currents in our analysis and work with forward 
transmission paths only. 

In the reconstruction of the sound speed perturbation field, the correlation 
covariance function of the field (or its discretized version, i.e., the covariance 
matrix Cs.) needs to be specified. We assume that the sound speed perturbation 
field is homogeneous and has a gaussian shape so that we can specify the 


correlation between the field at two different points in an analytical form as 











2 [fey fey] 
Cov(Ax,Ay,Az) = 65,€ ee Bre), 


where Ax is the horizontal separation between the two points in the x-direction, 
Ay is the horizontal separation in the y-direction, and Az is the vertical 
separation. The correlation length, Lx, Ly, and Lz, determine the correlation 
scale of the field. TABLE 4-1 summarizes the various sets of correlation lengths 
(Lx, Ly, and Lz) that we used to simulate the sound speed perturbation fields of 


different scales. For all the simulated fields, we use Osc = 5 m/s. 
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Figure 4-1: The CTD Data of Mooring #1 From the Deployment Cruise 
(Worcester and Howe, 1989). 
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TABLE 4-1 : THE SIMULATED OCEAN VOLUME. 


Ocean volume L, (km) L, (km) L, (km) 
Eddy204 20 20 0.4 
Eddy304 30 30 0.4 
Eddy404 40 40 0.4 


In order to perform simulation inversions in the computer, we need to 
generate a Set of travel time data for input to the estimator. A normal random 
deviate generator is used to generate realizations of a random process having a 
specified covariance. This random deviate generator is used to simulate the sound 
speed perturbation and noise fields. The simulated fields are then combined using 
the model (given by Eq. (2.9)) to give the simulated travel time data. A block 


diagram of the process is shown in Figure 4-2. 


Random 
Coc deviator 6c — | A ————» oT 
generator 


Random 
deviator a 
generator 





1o 
Ic 


Figure 4-2 : The Block Diagram of Travel Time Generation. 


The generated travel time data with additive random noise are the input to 


ue 
the estimator. An estimate of the sound speed field, Oc is the inverse result, 


which depend on the matrices C;,., C, and A in addition to the data themselves. 


The sensitivity of the system can thus be evaluated for changes in correlation 


a) 


length of the sound speed perturbation field, changes in noise level, changes in 
array geometry due to mooring failure, etc. To simulate noise-free 
measurements, we uSe a very Small noise level with an rms value of o, = 0.1 ms. 

The optimal estimate is obtained when the correlation length used for 
inversion is exactly the same as that actually present in the ocean volume (that is, 
when the a priori covariance Cs, is correct). In fact, because the a priori 
covariance is never exactly known, optimal estimates are generally hard to 
obtain. The quality of suboptimal estimates can be evaluated by studying the 
effect of correlation length uncertainties AL,, AL,, and AL, on the performance 


of the estimator. 


B. SYSTEM PERFORMANCE 
1. RMS Error Analysis 

The local RMS error as a function of box index j in the estimate (that 
is, the RMS error at the jth box) can be obtained by evaluating the square root of 
the jth element along the diagonal of the estimate error covariance matrix C,. 
This local RMS error gives a picture of how the errors are distributed spatially. 
The square root of the spatial average of the mean square errors at each of the 
boxes, Of, was calculated to get an idea of how much the local error varies over 
the whole ocean volume being studied. 

Figure 4-3 shows the contour plots of the local RMS estimate error for 
one set of system parameters in each of the five vertical layers. The RMS error 
for this case 1s approximate from 1 m/s to 2.5 m/s (or 20% to 50% compared to 
the 5 m/s signal level) and shows gradual dependence on horizontal or vertical 
position inside the perimeter of the array. The error outside the array's 


perimeter increases rapidly as distance from the array increases. Obviously, we 
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Figure 4-3: RMS Error (m/s) at Each Layer. The contour level is 0.5 m/s. 
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cannot obtain accurate estimate of the sound speed perturbations for points 
outside the array. There is little difference between the RMS errors in the 
individual layers because each vertical layer are well sampled by the eigenrays 
whose turning points occupy every layer. Because of the rather weak vertical 
dependence, in the following discussion we present only the analysis for the first 
layer. 

The local RMS error maps for the first layer for different horizontal 
correlation lengths are shown in Figure 4-4. The system parameters are identical 
to those of Figure 4-3, except that the horizontal correlation length is varied 
from 20 km to 60 km in a 10 km step. We see that a wider covariance matrix 
(i.e., a longer correlation length) results in a lower estimate error. This 
observation is not surprising, since a longer correlation length means that there 
is more gradual variation in sound speed perturbation with position, and that 
neighboring boxes are more correlated. An increased correlation reduces the 
number of degrees of freedom in the solution, and thus giving a better 
determined solution. 

In Figure 4-5, the maps are generated using the system parameters 
identical to those used to generate Figure 4-3, except that the vertical correlation 
length is varied from 0.2 km to 0.6 km in a 0.1 km step at a fixed horizontal 
correlation length of 40 km. There is less change in RMS error over the entire 
range of correlation lengths for this case than for the case shown in Figure 4-4. 
This is because the eigenrays sample the vertical layer more adequately than the 
horizontal sections. 

The effect of an increased (or decreased) noise level on the estimate are 


Shown in Figure 4-6. A higher noise level results only in a slightly higher RMS 
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Figure 4-4: Estimate RMS Error (m/s). The contour level is 0.5 m/s. The 
vertical correlation length L, is fixed at 0.4 km. 


4] 


(a) Tomography Array (d) L,=0.4km 
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Figure 4-5: Estimate RMS Error (m/s). The contour level is 0.5 m/s. The 
horizontal correlation length L, and L, is fixed at 40 km. 
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(a) Noise Level = 0.1 ms 
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Figure 4-6: Estimate RMS Error at Different Noise Levels. The contour 
level is 0.5 m/s. 
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error. The estimate is thus not as strongly affected by changes in noise as it is by 
changes in the ocean correlation length. 

Figure 4-7 shows the effect of array element failure on the RMS error 
(the array configuration is superimposed on the contour plots to give an idea of 
how the estimate depends on the array configuration). Within the perimeter of 
the remaining "good" elements, the RMS errors are only about 25% (0.5 m/s) 
higher than those given by the full array (2.0 m/s). However, the errors outside 
the area covered by the good elements increase rapidly to 100%. The failure of 
array elements has a very pronounced effect on the estimate, especially when two 
or more elements fail. The system essentially loses its ability to give accurate 
estimates in areas that do not have rays passing through them. 

A single measure of RMS error can be evaluated by calculating a spatial 
“average of the individual local errors. We use Og for such a global measure 
and it is computed by assuming the individual mean-square errors at all the box 


locations and then taking the square root of that sum. Figures 4-8 and 4-9 show 


Of as a function of the ocean horizontal and vertical correlation lengths, 
respectively, for various array geometries and noise levels. In all cases, Og was 
found to decrease with increasing L, and Ly (see Figure 4-8). However, O¢ is not 
as sensitive to changes in L, as it is to changes in L, and L, (see Figure 4-9). This 
is due to the fact that the vertical structure is sampled adequately whereas the 
horizontal structure is not. In view of this, in the following discussion, we will 
only discuss the sensitivity of the results of our resolution analysis with the 


vertical correlation length fixed at 0.4 km. 
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Figure 4-7: Estimate RMS Error (m/s) in Mooring Failure Cases. The 
contour level is 0.5 m/s. Ly= L, = 40 km, and L, = 0.4 km. 
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Figure 4-8: Square Root of Spatial Average Mean Square Error vs. 
Horizontal Correlation Length. 
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Figure 4-9: Square Root of Spatial Average Mean Square Error vs. Vertical 
Correlation Length. 
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2. Resolution Analysis 

Two measures of system resolution were defined in the last chapter. 
One measure is the minimum resolution length, which gives the local resolution 
at each box. The minimum resolution length is essentially the size of the smallest 
ocean feature which can be resolved by the array. A large minimum resolution 
length indicates a poor resolving power. Figure 4-10 shows the minimum 
resolution length #, in the x direction as a function of L, and L,. We see that 
for L, and L, greater than 30 km, #, 1s relatively constant over the interior of 
the array and in most of the region of interest. On the other hand, if L, and Ly 
are less than or equal to 30 km, # varies heavily on both x and y and becomes 
very large in regions containing no y-oriented ray paths. The behavior of 4, the 
minimum resolution length in the y direction, is analogous to the behavior of H, 
and is shown in Figure 4-11. In general, the average minimum resolution length 
is approximately 30 km inside the monitored region as shown in Figure 4-10 and 
Figure 4-11. 

Figure 4-12 shows the behavior of +, and #, as noise level change 
while the horizontal and vertical correlation lengths are fixed at 40 km and 0.4 
km, respectively. We see that varying the noise level has little effect on #, and 
H,. Figure 4-13 shows how the failure of various array elements affects H,. A 
diagram of the array configuration is superimposed on each plot in Figure 4-13 
to show the connection between array geometry and H,. Despite the element 
failures, H, remains relatively constant inside the area covered by the remaining 
array elements, but it increases rapidly as we move outward from the perimeter. 
Although there is some resolving power outside the array perimeter, it is highly 


inadequate. 
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Figure 4-10: Horizontal Minimum Resolution Length H,. The contour level 
is 10 km. L, Its fixed at 0.4 km. 
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(a) Tomography Array (d) L,, Ly = 40km 
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Figure 4-11: Horizontal Minimum Resolution Length H,. The contour level 
is 10 km. L, is fixed at 0.4 km. 
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Figure 4-12: Horizontal Minimum Resolution Lengths at different Noise 
levels. The contour level is 10 km. L, = Ly = 40 km, and L, = 0.4 km. 
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Figure 4-13: Horizontal Minimum Resolution Lengths in various Mooring 

Failure cases. The contour level is 10 km. L, = Ly = 40 km, 

and L, = 0.4 km. 
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The second measure of system resolution is the resolution spread (Eq. 
3.16), which is the square of the Frobenius norm of the difference between the 
resolution matrix and the identity matrix. The resolution spread is thus a single 
scalar quantity which describes the resolution of an estimator over an entire 
region. Figure 4-14 shows resolution spread as a function of horizontal 
correlation length for each of several noise levels and array geometries. We see 
that the resolution spread does not significantly depend on changes in correlation 
length, except when noise level 1s increased. That is, system resolution is a lot 
more sensitive to changes in correlation length when the system 1s not noise free. 
Generally, the resolution spread depends most strongly on the array geometry at 


any fixed correlation length. 
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Figure 4-14: Resolution Spread vs. Horizontal Correlation Length. 
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C. SENSITIVITY OF INVERSE SOLUTION 

Synthetic travel time perturbation data OT (refer to Figure 4-2 for the 
generation of ST) were used to reconstruct the sound speed perturbation fields 6c 
simulated in the computer. The sensitivity of the estimate to uncertainty in the 
correlation length specified for inversion is examined in this section by 
comparing the suboptimal estimates to the optimal estimates. The optimal 
estimates are those derived using an exact covariance of the field whereas the 
suboptimal ones are consequences of inexact covariances. 

Figure 4-15 shows the simulated field Eddy404 (refer to TABLE 4-1 for the 
simulation parameters), its optimal estimate, the associated RMS error as well as 
the difference between the estimate and the simulated field. The estimate error 
&e-5C is low (from 0 to +2 m/s, which is from 0% to 40% compared to a signal 
level of 5 m/s) over most of the area inside the array. The error is larger in the 
left edge and corners. Figures 4-16 compares the estimate errors of some 
suboptimal estimates generated using correlation lengths differ from the true 
one. The difference between the assumed horizontal correlation lengths for 
inversion and the true lengths are denoted by AL, and AL, in the figure. 
Obviously, the effect of a positive correlation length uncertainty seems less 
harmful than that of a negative one when the actual correlation length is 40 km. 

Figure 4-17 shows the effect of uncertainty in the noise level (i.e., in the 
noise covariance matrix C,) on the estimate. The maps on the left are the 
estimates of the Eddy404 field with various noise level uncertainties Ao,. On the 
right the associated errors are displayed. Generally, a higher noise level 
uncertainty gives a higher estimate error, although the differences in the errors 


are quite minimal. 
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Figure 4-15: A simulated field, its optimal estimate, estimate RMS error, 
and difference between the optimal estimate and the simulated field. Units 


are in m/s. L,= Ly = 40 km, L, = 0.4 km, and AL,= AL, =AL,= 0. 
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Figure 4-16: Simulated Eddy404 field, suboptimal estimates, and the 
corresponding true error fields. The contour level is 4 m/s and 2 m/s for 


S¢ and Sc — Sc, respectively. L, = Ly = 40 km, and L, = 0.4 km. 
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Figure 4-16: (Continued from last page). 
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Figure 4-17: Estimation (m/s) with Noise Level Uncertainty. The contour 


level is 4 m/s and 2 m/s for $¢ and 5c ~ dc, respectively. Ly = Ly = 
40 km, and L, = 0.4 km. 
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To gain insight into the global performance of the estimator, we calculated 


O;, the square root of the spatial average of the local squared errors. Figure 4-18 
shows 6; as a function of horizontal correlation uncertainty for each of the three 


simulated ocean volumes having horizontal correlation length of 20, 30 and 40 


km, respectively. 
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Figure 4-18: Square Root of Spatial Average Square Error vs. 
Horizontal Correlation Length Uncertainty. 


It is interesting to note that, for assumed horizontal correlation lengths of L, 
and L, both equal to 1 km’, the error is 3.4 m/s for all the three cases estimated. 


Furthermore, Ox is less sensitive to a positive correlation length uncertainty than 


* Correlation length of zero implies no a@ priori information. Our software 
does not permit the use of zero correlation length, so we used correlation length 
of 1 km to approximate the case of no a priori information. 
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to a negative correlation length uncertainty when the actual ocean correlation 
length is long (e.g. longer than 20 km, see the curves associated with ocean fields 
Eddy304 and Eddy404), while it is more sensitive to a positive uncertainty than 
to a negative uncertainty when the actual correlation length is short (e.g. 20 km, 
see the curve associate with Eddy204). 

In an ocean with small eddies, Figure 4-18 suggests that using a correlation 
length longer than the actual one could result in a large increase of estimate 
error. Whether or not the estimate is sensitive to an inexact and longer 
correlation length used really depends on the resolution of the system. We see 
from Figure 4-10 that the minimum horizontal resolution length of the 
Greenland Sea array 1s approximately 30 km, almost regardless of what the 
correlation length of the field is. This implies that the array is unable to resolve 
ocean features smaller than 30 km. In other words, an ocean field with a 
correlation length shorter than 30 km (e.g. 20 km) can not be monitored 
adequately. In that case any additional a priori information, even though wrong, 
is totally absorbed by the estimator for solution construction, thus leading to an 
uncertainty sensitive estimate. In contrast, for an ocean field with a correlation 
length longer than 30 km, the array turns into an adequate system. The 
corresponding estimator now has the ability to reject extraneous information. 


The result is an uncertainty insensitive estimate. 


TABLE 4-2 shows the value of o¢ for each suboptimal estimate, as well as 


the percent difference in 6g between each suboptimal estimate and the optimal 


estimate for the ocean field Eddy404. By selecting a value for the maximum 


acceptable percent difference in Og between the suboptimal estimators and the 


optimal estimator, we can compute from the curves in Figure 4-18 the maximum 


a), 


allowable uncertainty in correlation length. TABLE 4-3 shows, for each 

simulated field, the allowable range of uncertainty in correlation length if the 

difference between 6; for the suboptimal and optimal estimators is not to exceed 

ten percent. 

TABLE 4-2 : THE ESTIMATE ERROR AND @% DIFFERENCE TO 
THE OPTIMAL ESTIMATE FOR EDDY404. 


AL,, AL, (km) Oz % difference 
-30 3.30 (Rees, 
-20 2.50 32S 
-10 2.06 D511 
10 1.88 0.00 
20 Lee 9.63 
30 2.10 11.70 


TABLE 4-3 : TEN % DIFFERENCE ALLOWABLE UNCERTAINTY 


RANGE. 
L,, L, (km) Allowable AL, and AL, Range (km) 
20 km -20.0 ~ 15.6 
30 km -11.8 ~ 23.7 
40 km -10.4 ~ 25.8 


TABLE 4-3 quantitatively confirms our earlier observation from Figure 4 
18 that the estimate error is less sensitive to a positive uncertainty in correlation 
length than to a negative one in ocean volumes containing large structures, and 


the result is reversed for an ocean volume containing small structures. 
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V. CONCLUSION 


A. ESTIMATOR PERFORMANCE 

In order to obtain a unique estimate of the sound speed perturbation field, a 
priori information in the form of a sound speed perturbation covariance matrix 
is used in the Gauss-Markoff estimator. As discussed in Chapter 4, given that the 
sound speed perturbation field has a gaussian shape correlation, the optimal 
estimate for the field is obtained when the assumed correlation length (i.e., the 
correlation length used to calculate the input covariance matrix C;.) is equal to 
the actual correlation length present in the ocean. A primary goal of this thesis 
is to evaluate the effect on our estimator when the assumed correlation length is 
not equal to the actual one. As this happens, the estimate becomes suboptimal. 
When the actual correlation length is not exactly known, it has been suggested by 
Cornuelle (1985) and Chiu (1987) to use a conservative assumption (J.e., a small 
correlation length) so as not to "assume too much" about the sound speed 
perturbation field. 

From our simulation study we found that the optimal estimate for the sound 
speed perturbation field typically has an RMS error between | and 2 m/s (i.e., 
20% to 40% comparing to the 5 m/s signal level). For the suboptimal case, the 
estimate is actually less sensitive to a positive correlation length uncertainty (1.e., 
the assumed correlation length is longer than the actual one) than to a negative 
uncertainty when the actual correlation length is longer than 30 km. On the 
other hand, when the actual correlation length is shorter than 30 km (for 


example 20 km) the estimate becomes more sensitive to a positive uncertainty 
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than to a negative uncertainty. The reason for this switching of behavior at a 
correlation length of 30 km has something to do with the resolution of the 
Greenland Sea array. In our resolution analysis we found that the Greenland Sea 
array has a resolution length of about 30 km. Therefore, for an ocean field with 
a correlation length small than 30 km (for example 20 km), the use of a 
correlation length longer than the actual one for inversion would ingest 
extraneous information into an information hungry estimator. This estimator 
basically accepts all the wrong information. It is thus preferable to be 
conservative and use a correlation length which is likely to have a negative 
uncertainty when the system resolution is inadequate for measuring the expected 
scale. On the other hands, if we know from resolution analysis that the resolution 
is adequate for a particular ocean region, a positive correlation length 
uncertainty is acceptable in this case. 

By specifying the maximum acceptable percent difference in Og (which is 
essentially a spatial average of local estimate errors) between the suboptimal and 
optimal estimators, we arrive at one possible criterion for designing the 
estimator. If this difference, in the case of an expected ocean correlation length 
of about 40 km, is required to be less than ten percent (for example), the 
estimator can tolerate any correlation length uncertainty between -10.4 and 25.8 
km, a spread of 35.6 km. 

We have also shown that a higher noise level results in basically unchanged 
RMS error and resolution. This result suggests that the estimate error is not 
dominated by the random error, but rather by the bias error arising from the 


fact that the tomography problem is underdetermined. This insensitivity to 
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random noise is one benefit of using a Gauss-Markoff estimator, which always 
tries to minimize the effect of random noise (Chiu et al, 1987). 

The failure of array elements has a very pronounced effect on the RMS 
error and resolution, especially if two elements fail. The comparison of cases 
simulating the failure of one or more elements indicates that the RMS error 
becomes very large and resolution becomes very poor in areas that no longer 
have rays passing through them. However, in regions still containing acoustic 


rays, the RMS error is only 25% higher than that of the full array case. 


B. RECOMMENDATIONS FOR FUTURE IMPROVEMENT 

All numerical simulations were performed on a DEC MicroVAX. Each run 
required roughly five and one-half hours of CPU time. A machine of greater 
computational power and larger memory size would allow us to divide the ocean 
volume into finer spatial meshes for improved analysis. 

We can derive statistical information concerning the vertical ocean structure 
from historical data using the empirical orthogonal function (EOF) approach of 
Cornuelle (1983, pp. 139). EOF analysis has been widely applied in research 
fields other than tomography. In EOF analysis, the vertical structure of the ocean 
is represented by a set of orthogonal vectors. These vectors can be derived from 
a singular value decomposition of a matrix containing historical data from 
hydrographic surveys. The EOF method gives a priori information about the 
vertical structure which is more realistic than that provided by the Gaussian 
shape correlation used in this study, and is recommended for use when the actual 


tomographic data from the Greenland Sea Project become available. 
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